home *** CD-ROM | disk | FTP | other *** search
/ Resource for Source: C/C++ / Resource for Source - C-C++.iso / codelib8 / v_10_05 / 1005032a < prev    next >
Encoding:
Text File  |  1995-11-01  |  2.1 KB  |  52 lines
  1. #include <math.h>
  2.  
  3. /* Adaptive quadrature routine.
  4.    Called with:
  5.         lft:            left endpoint
  6.         rt:             right endpoint
  7.         f:              pointer to function to integrate
  8.         err:            pointer to allowable error
  9.  
  10.    Returns
  11.         return value    evaluated integral
  12.         err             pointer to upper limit of actual error
  13. */
  14.  
  15. double adaptive (double lft, double rt, double (*f)(double x), double *err)
  16. {
  17.    double h;                           /* Width of this interval. */
  18.    double p1, p2, p3, p4, p5;          /* Value of function at the 5 points. */
  19.    double error;                       /* Actual error on this interval. */
  20.    double err1, err2;                  /* Errors of the sub-intervals. */
  21.    double int3,                        /* Value of the integral with Simpson */
  22.           int5;                        /* 3- and 5-point approximations. */
  23.  
  24.    h = rt - lft;                       /* Get the width of the interval. */
  25.                               /* Get the 2 approximations to the integral. */
  26.  
  27.    p1 = (*f)(lft);
  28.    p2 = (*f)(lft + h/4.0);
  29.    p3 = (*f)(lft + h/2.0);
  30.    p4 = (*f)(rt - h/4.0);
  31.    p5 = (*f)(rt);
  32.  
  33.    int3 = h * (p1 + 4.0 * p3 + p5) / 6.0;
  34.    int5 = h * (p1 + 4.0 * p2 + 2.0 * p3 + 4.0 * p4 + p5) / 12.0;
  35.  
  36.    error = fabs (int3 - int5) / 15.0;     /* Find the actual error. */
  37.    if (error < *err)                      /* If within tolerance, */
  38.    {
  39.       *err = error;                       /* copy the actual error to caller, */
  40.       return (int5);                      /* return the more accurate approx. */
  41.    }
  42.    else                                   /* Otherwise, */
  43.    {
  44.       err1 = err2 = *err / 2.0;           /* cut the allowable error in half, */
  45.       int5 = adaptive (lft, lft + h/2.0, f, &err1);     /* integrate the two */
  46.       int5 += adaptive (lft + h/2.0, rt, f, &err2);     /* half-regions, */
  47.       *err = err1 + err2;                 /* sum the real errors, */
  48.       return (int5);                      /* and return the integral. */
  49.    }
  50. }
  51.  
  52.